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Abstract 

Reversible bipolar nano-switches that can be set and read electronically in a solid-state two- 
terminal device are very promising for applications. We have performed molecular-dynamics simu- 
lations that mimic systems with oxygen vacancies interacting via realistic potentials and driven by 
an external bias voltage. The competing short- and long-range interactions among charged mobile 
vacancies lead to density fluctuations and short-range ordering, while illustrating some aspects of 
observed experimental behavior, such as memristor polarity inversion. 
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I. INTRODUCTION 



Nanoscale metal /oxide /metal two-terminal non-linear circuit elements and switches 
would be extremely attractive for dense memory, logic, neurocomputing etc if they could 
scale well in size, power, driving voltage, and can perform without much fatigue in a repeat- 
able and uniform (across e.g. the memory matrix) manner. Some of these issues are not 
well understood although the hysteretic behavior of such materials, especially thin films of 
Transition Metal Oxides (TMO) like Ta20s, Nh^Os, TiC>2, NiO, CU2O, and Group III and 



IV oxides (AI0O3, SiO x ), in metal-insulator-metal (MIM) vertica 



111 ], polymer 



12 



devices have been studied 



13| . and correlated oxide 



for decades Il|j8||. The diverse semiconductor! 9 

systems lol. Il4j-|l7( exhibited switching under electric pulsing. In recent years, interest in 
various oxide-based systems switchable by an electric pulse has grown quite dramatically 
leading to important breakthroughs and exposing various fundamental problems related to 
mechanisms of bipolar (driven by alternating positive and negative biasing to close the I-V 
hysteresis loop) and unipolar (driven by the bias of one polarity) switching behavior 18H20]. 
In particular, it was realized that bipolar switches are analogous to the 'memristor', a fourth 
passive circuit element originally postulated by Leon Chua in 1971 [18]. 

There are challenges in understanding and controlling the coupled electronic and ionic 
kinetic phenomena that dominate the behavior of oxide switching devices like Pt/TiO^/Pt, 
which is an exemplar memristor (resistor with memory) 19] . It has been demonstrated un- 
ambiguously that bipolar switching involves changes to the electronic barrier at the Pt/Ti02 
interface due to the drift of positively charged oxygen vacancies under an applied electric 



field 



191 ] . Various direct measurements revealed formation of localized conducting channels 



in TiO"2: pressure modulated conductance microscopy 



21 



221 ]. conducting atomic force 



microscopy ( AFM) [23!] , scan ning transmission x-ray microscopy 23, |25) , and in-situ trans- 



mission electron microscopy 



261 ] . On the basis of these measurements, it became quite 



clear that the vacancy drift towards the interface creates conducting channels that shunt, 



or short-circuit, the electronic barrier and switch the device ONjl9j. The drift of vacancies 
away from the interface breaks such channels, recovering the electronic barrier to switch the 
junction OFF. More recently, Kwon et al. 26( have performed the direct cross-sectional 
TEM studies of the unipolar resistive switching of Ti02, revealing the presence of nanoscale 
Magneli phase Ti 4 07 conductive channels following device turning on and imaging the fila- 
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ment region in TEM. A Ti407 phase was confirmed from the temperature dependence of the 



25j investigated the bipolar mode of TiC>2 



conductance. Concurrently, Strachan et al. [24, 
switching, using non-destructive spatially-resolved x-ray absorption and electron diffraction 
that allows nanometer scale studies of the associated chemical and structural changes of 
a functioning memristive device. They have observed that electroforming of the device is 
accompanied by forming an ordered Magneli phase channel within the initially deposited 
(amorphous) Ti02 matrix. 

The observation of a Magneli crystallite in a titanium dioxide matrix shows that electro- 
forming Pt/Ti02/Pt memristive system is related to a localized partial reduction of titanium 
dioxide and crystallization of a metallic conducting channel Ti 4 07. Inside a titanium diox- 
ide matrix, the Magneli phases are thermodynamically favored over a high concentration 
of randomly distributed vacancies, and thus they can act as a source or sink of vacancies 
in the matrix material depending on the electrochemical potential within the device 271]. 
The application of an electrical bias can control vacancy motion in and out of this sub-oxide 
phase, modulating a transport barrier and leading to the dramatic conductivity change. It is 



worth noting that hysteretic switching in some perovskite oxides monitored by TEM did not 



show any indication of conducting channel formation 



Therefore, two types of models 



are usually considered (i) in-plane homogeneous and (ii) in-plane inhomogeneous changes 
accompanying switching in the material 7j. We shall study the kinetics of oxygen vacancies 
in a generic situation like TiC>2 by way of molecular dynamics to visualize processes taking 
place at the atomic scale that are reminiscent of actual device behavior and whose origin is 
in competing vacancy-vacancy interactions. 

Another important aspect is that the forming step for the channel formation in systems 
like TiC>2, NiO, VO2 is accompanied by local heating that is witnessed by the local emergence 
of high-temperature phases and observed by thermal microscopy. It is worth noting that 
bipolar switching is apparently field-driven and is very different from unipolar switching 
that is power-driven. As a result, the bipolar switches may consume much less power 
than the unipolar switches. At the same time, local heating does accompany and seems 
to be important in both processes. Thus, fresh TiC>2 samples have an amorphous layer of 
titanium dioxide. To make the device switch in a repeatable manner, a forming step of 
rather large voltage pulse is usually required to create a vacancy-rich region (forming is 
not required if e.g. a vacancy- rich layer is provided intentionally by fabrication). During 
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this process, the Ti02 anatase phase forms near the conducting channel, which nominally 



requires temperatures over 350° C 



251 ] . Studying local thermal effects during switching 



in Ti02 provides strong evidence for local heating 29j. Unambiguous evidence for local 



leating accompanying conducting channel formation comes from third harmonic generation 



30| ; low resistance states (LRSs) in NiO showing unipolar switching were strongly nonlinear 
with variations in the resistance R as large as 60%, which was most likely caused by the 
Joule heating of conducting channels inside the films. One can note at this point that since 
the unipolar switching is power driven, therefore, it may be problematic to use it because 
of power requirements. 

Given the above, the microscopic understanding of the atomic-scale mechanism and iden- 
tification of the material changes within an insulating barrier appears to be invaluable for 
controlling and improving the memristor performance. Within the volume 10x10x2 nm 3 of 
perspective nano-memristors the number of oxygen vacancies could be as small as a thou- 
sand, so that the conventional statistical approach for dealing with many-particle systems 
might fail. Many important phenomena such as dynamic phase transitions, as well as compe- 
tition between thermal fluctuations and particle-particle interactions in stochastic transport, 
can be investigated using the Molecular Dynamics (MD) simulations of the Langevin equa- 



tions describing the thermal diffusion and drift of individual interacting particles 31J. Here 
we present some preliminary results of such simulations for a model memristor. 

II. MODELING A MEMRISTOR 

In a semiconducting titanium dioxide, oxygen vacancies play an important role in the 
electrical conduction. A stoichiometric Ti02 crystal is a colorless and transparent insulator. 
However, heated in vacuum, it loses a part of its oxygen, or is reduced, and becomes dark blue 
in color and electrically conductive. At high temperatures about 1000 C ionic conduction 



becomes observable originating f rom oxygen and tit« vacancies U while at ambient 
and lower temperatures small polarons carry the current [33!]. The reduced rutile shows an 
anomalous electrical property under a DC electric field of several tens of V/cm [32[. When 
a DC electric field was applied to the bulk sample of reduced Ti02, a single pulse of electric 
current was observed together with a constant background current. This pulsed current has 
been attributed to migration of oxygen vacancies in reduced Ti02 under a DC electric field. 
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FIG. 1: (Color online) Model oxide memristor with a reduced oxygen layer near one of the metallic 
leads (the length scale is arbitrary) 



The oxygen migration is visually observable under a DC electric field. Namely, when a DC 
electric field was applied to the sample after repeated change of its polarity, the color of the 
sample around the negative electrode changed into dark blue due to the higher density of 
oxygen vacancies, while the region around the positive electrode lost color 32J. In order to 
explain these features, Miyaoka et al. 32] have proposed that the oxygen vacancies form 
clusters under a voltage, which may create crystalline phases different from rutile, namely 
Magneli phases, as recently observed in nanoscale probes 25l. \26\. 



32| and 



Based on the observations of oxygen vacancy migration and clustering in bulk 
nanoscale 125], 



26| samples of Ti02 induced by an electric-field, we model a memristor as 



shown in Fig. HJ In our model, there is a reduced rutile thin layer Ti02_ x near one of the 
metallic electrodes stabilized by the Coulomb mirror potential. Vacancies from this layer 
can drift toward the opposite electrode in an electric field. 

The simplest approach to the vacancy diffusion is using a drift-diffusion equation for the 
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density of vacancies n(x,t), 

~d 2 n(x,t) F(t)dn{x,t) 



dn(x,t) = 
dt 



dx 2 ksT dx 



(1) 



which neglects the particle-particle interaction. Here D is the diffusion coefficient, T is the 
temperature, and F{t) is the electric field. The diffusion equation (TTJ can be mapped onto 
a simple Langevin equation 

F*{t) + VDZ$(t), (2) 



dt 

where a = x, y, z Cartesian components, £ is a stochastic force of zero mean and 5-correlated 
in time, and we accounted for anisotropic diffusion (in Ti02 rutile diffusion along the c-axis 
is much faster than in ab-plane, for instance). 



With the use of the standard Euler method 



3l| , we have simulated one thousand vacancies 



(i = 1, 1000). The result is shown in Figf5]for a rectangular electric pulse, F(t). There is 
no clustering in this simple diffusion approximation, but rather uniform propagation of the 
reduced layer accompanied by its widening in the propagation direction, which is the same 
as the analytical result in this case. 

III. MOLECULAR DYNAMICS OF INTERACTING VACANCIES 

Now let us simulate the stochastic transport of vacancies fully taking into account their 
interaction and boundary conditions. We consider an open system of N point-like Brownian 
vacancies interacting with each other via the pairwise potential W, with a possible inclusion 
of the substrate through some periodic and random potentials U, and with internal and 
external fields corresponding to the time-dependent deterministic force F. The environment 
exerts an independent Gaussian random force, £ on each particle with zero mean and inten- 
sity controlled by the temperature. In the overdamped regime (where inertia is negligible 
compared to the viscous damping), the Langevin equation describing the drift-diffusion of 



the i-th particle is [31] : 



, = k(*u t) - E awi Z Xs) + v^WI» ( 3 ) 



The unit of £j is s 1 ^ 2 . The physical unit of length is the thickness of the memristor, L, 
(about or less than 5 nm). Then the physical unit of time is L 2 /D, where D = hsT/i]. 
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t=0.01 



t=0.05 



t=0.07 



FIG. 2: (Color online) Simulated vacancy motion under the electric pulse governed by the simplest 
diffusion equation ([TJ showing no clustering (length and time scales are arbitrary) 

(we take the largest D for anisotropic TiC^). The diffusion coefficient D should be taken at 
a local temperature, but we have assumed some net temperature across the active area in 
order to simplify matters in the present simulation. Using the dimensionless time r = tD/L 2 
and dimensionless coordinates Tj = Xi/L, the equations ([3]) can be written as 



where V(ri,r) is the electric-field potential, and the components of the (dimensionless) 
random force, £f (a = x,y,z), satisfies the fluctuation - dissipation relation (£f(0)o( T )) = 
8{r)8 a fi5ij. Diffusion in Ti02 is strongly anisotropic, thus the diffusion along the c-axis 
and in-plane will be described by two equations with different diffusion coefficients in more 
accurate simulations. 

The set of equations (jlj) is equivalent to an infinite BBGKY (Bogoliubov - Born - Green - 
Kirkwood - Yvon hierarchy) set of equations for multi-particle distribution functions. With 
no particle-particle interaction it reduces to the drift-diffusion equation ([TJ). Also, in a mean 




g[U(r t ,r) + [/(r t )]/(A; J3 T) 
drf 
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field approximation, one obtains an additional drift term in the diffusion equation, which 
however does not lead to clustering. 

The O 2- vacancy-vacancy interaction potential, which is crucial for their clustering and 
phase transformations, can be modeled as 

W(R) = Aexp(-R/a) - B/R e + e 2 /(ire eR), (5) 

where the short-range repulsive and attractive parts are represented with A = 22764.0 eV, 



a = 0.01490 nm, and B = 27.88 x 10 _e eV nm [34|, and the long-range Coulomb repulsion 
is the last term. More generally the interaction parameters A, B, a vary from one oxide to 
another, and the dielectric constant e may also vary from sample to sample. 



To estimate the diffusion coefficient one can use an empirical equation of Ref. 



35|, 



D [cm 2 /s] = (1.03 x 10" 3 ) exp (-Ejk B T) , (6) 

where the activation energy for oxygen vacancy diffusion has a poorly known value, with 
estimates ranging from 0.5 eV to about 1.1 eV. The lower boundary, E v = 0.5eV, corresponds 
to D pa 3.44 x 10~ 12 cm 2 /s at the room temperature. For higher activation energies, the 
vacancies would be immobile, as they are known to be at room temperature. We do know, 
however, that they move during switching. In this regard, it is worth noting that (i) another 
important parameter is the local temperature during switching and (ii) strong local electric 
field concentration at the tips of conducting channels/filaments. Some estimates indicate 
that the local temperature may reach values like 650K, in which case the diffusion coefficient 
would be much larger, obviously. The local field concentration near tips of conducting 
channels/filaments is another variable that will very strongly facilitate the growth of the 
channel and needs to be simulated self-consistently along with local heating. 

Given the uncertainties with ionic diffusion in memristors, at this stage we would restrict 
ourselves to qualitative and semi-quantitative analysis. The following parameters are chosen 
at this stage more for computational convenience than to model actual devices. In actual 
devices the voltage drop across the active region may be moderate, on the order of IV, 
and the pulse duration for switching may be just a few microseconds. With L = 5 nm 
and such diffusion coefficient the time unit is L 2 /D 0.7 ms. The electric pulse rising 
time 100 ns would correspond to r r = 1.4 x 10~ 4 and an pulse length 10 /xs corresponds to 
T im P = 1-4 x 10~ 3 . The pulse amplitude of 5 V would then yield the dimensionless force 




FIG. 3: (Color online) Two successive runs (upper and lower panels, respectively) of the toy-model 
simulations with periodic BC and parameters described in the text. The bottom panel illustrates 
the "polarity inversion" of the vacancies from bottom to top electrode during successive pulsing. 



Analogous behavior has been observed experimentally in Ref. 



36]. 



F = eV/(ksT) w 200 at room temperature. The characteristic particle speed is v — eV/Lrj, 
so that the dimensionless collision time is t co i = a/LF « 1.5 x 10~ 5 , and the number 
of computing steps is l/r co ; 70000, representing a challenge for MD simulations. The 
characteristic collision distance R c , where the collision force changes its sign from attraction 
to repulsion, is in fact several times larger than a, which makes the number of MD steps 
manageable. Also, decreasing the electric pulse amplitude makes the number of the MD 
steps smaller. We will present simulations for more realistic parameters elsewhere, with the 
inclusion of local thermal effects directly into simulations. 
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FIG. 4: (Color online) The same as in Fig. [3] with a trend towards formation of a percolation path 
along the chains of oxygen vacancies. 

IV. MD SIMULATION OF A TOY MEMRISTOR: CLUSTERING OF VACAN- 
CIES 

In view of those uncertainties, we present here the MD simulations of a toy memristor 
with relatively small number of vacancies, N = 75 (Fig. [3]) and 150 (Fig. |3J), leaving 
simulations of real devices for future studies. Initially, we have randomly placed all the 
vacancies near the bottom of a toy sample (see Fig. [3] for initial distribution panels) and 
then let these vacancies evolve according to Eq.(|4j) inside a rectangular box which mimics 
either a part (Figs. El Hj) or the entire (Fig. |5]) sample. We use the 2D simulation area 
with aspect ratio L y /L x = 2 (Figs. I3"f4"|) and 4/3 (Fig. [5]) and either periodic boundary 
conditions (BC) (Figj3l @| or impenetrable boundaries (Figj5j) along the x-direction. Note 
that periodic BC allow us to simulate a quite large sample using a rather small number of 
particles, while the impenetrable BC can elucidate the role of boundaries or ID defects in a 
real finite size sample. To simulate periodic BC, we include periodic images of vacancies with 
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FIG. 5: (Color online) Results of toy-model simulations of a system with impenetrable boundaries 
and parameters described in the text. As in Fig. [3l one can see clustering of Oxygen vacancies. 
Moreover, the sample boundary pins vacancies, favoring the formation of percolation paths at the 
edges of the sample. The left and right panels correspond to two successive runs with the same 
parameters but different initial vacancy distributions. The observed edge effect is very reminiscent 
of some memristors where switching areas tend to form near edges 25 ] . 



respect to vertical boundaries of the simulation box. We also incorporate opposite polarity 
charges by adding mirror images of vacancies with respect to the top and the bottom of the 
sample. Since we simulate a limited number of vacancies, one can refer to each particle in 
our simulations as a cluster of vacancies, where a cluster-cluster interaction described by the 
combination of the Lennard- Jones and Coulomb potentials acting with the force 

F(r) = 1 {l2E LJ [(r min /r) 12 - (r mm /r) 6 ] + E c r min /r) , (7) 

the relative strength of the Coulomb potential E c /Elj = 2 and the potential well is about 
30/c^T. This results in the position of the potential maxima r max m 2r min and the height of 
the potential barrier on the order of the depth of the potential well. The pulse strength is 

11 



taken to be about 10 times stronger than the maximum attracting force between vacancy 
clusters. 

The results of our toy-modeling look very promising. We see that pulsing in opposite 
directions initially pushes the vacancies towards the top electrode and then back, FigJSJ 
Pulsing twice in the forward direction moves the vacancies from the bottom to the top 
electrode, thus inverting the vacancy distribution in close analogy to what has been done 



experimentally 



361 ] . We observe a complex kinetics of the oxygen vacancies in all present 



simulations. In particular, we have observed fragmentation of vacancy distribution and 
formation of the short vacancy chains for both the periodic and impenetrable boundary 
conditions. One can envisage that in certain conditions those chains may form percolation 
paths (Fig. Hj), facilitated by e.g. pinning centers and/or boundaries of the sample. Further 
studies would elucidate the conditions under which a conducting channel may emerge in the 
sample driven by E-field, thus gaining more insight into the switching of memristors. 



V. CONCLUSIONS 



We have proposed a model for the kinetic behavior of oxide memristors and simulated its 
toy analog using the Molecular Dynamics Langevin equations. Our MD simulations reveal 
a significant departure of the vacancy distributions across the device from that expected 
within a standard drift-diffusion approximation, driven by interactions among the mobile 
ionic species (oxygen vacancies in Ti02, as a generic example). A realistic vacancy- vacancy 
interaction incorporating short and long-range repulsions along with a medium-range attrac- 
tion leads to clustering of vacancies into chains. This may shed new light onto conducting 
channel formation in memristors. We have found a significant effect of boundaries on the 
clustering patterns. MD experiments (before averaging over configurations) have the poten- 
tial to simulate conditions that are hard or even impossible to incorporate in the standard 
drift-diffusion models, such as the vacancy annihilation or generation (from the air), different 
shapes of boundaries, vacancy pinning (either periodic or random) and an inhomogeneous 
external field (due to edge effects etc.). One can also calculate not only the density distri- 
bution but binary distributions, velocity distributions etc., providing any desired statistical 
analysis of the distributions. Another important step would be to include thermal effects via 
local changes in the anisotropic diffusion coefficient. Competing interactions between mo- 
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bile ionic species are seen as a robust mechanism for vacancy concentration fluctuations and 
clustering, which shows that switching in real devices is a statistical or stochastic process. 
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